Mixed-Effects Models
Disturbance through human hunting activity can significantly impact prey species through both consumptive and nonconsumptive effects. The nonconsumptive effects of rabbit hunting on Northern Bobwhite (Colinus virginianus; hereafter, bobwhite) could cause an increased perceived risk of predation by bobwhite during rabbit hunting events may elicit anti-predator responses, such as reduced movement away from the safety of cover.
bobwhite <- read.csv('bobwhite3.csv')
bobwhite$ID <- as.factor(bobwhite$ID) #make ID a factor
p1<- ggplot(bobwhite, aes(x=HuntDay, y=HW_Dist, group=ID, color=ID, shape=ID)) +
geom_point(size=4, alpha=0.6, position = position_dodge2(width=.33, preserve = "total")) +
scale_y_continuous() +
#geom_line() +
geom_smooth(method = "lm", se = FALSE) +
labs(title="Risk Behavior in Bobwhite During Hunting Season", x= "Hunting Season Species", y = "Distance from Hardwood Forest Cover (meters)")+
theme_bw()+
scale_color_brewer(palette = "BrBG")
p1## `geom_smooth()` using formula 'y ~ x'
# ID: Unique ID given to each bobwhite covey tracked in chronological order.
# HuntDay: Denotes if it was a “Rabbit” or “Quail” (bobwhite) scheduled hunt day.
# HW_Dist: Distance in meters a bobwhite covey was from hardwood habitat.mixed_bobwhite <- lmer(HW_Dist~(HuntDay*ID)+(1|ID), data = bobwhite)## Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
## eigenvalue close to zero: 5.5e-09anova(mixed_bobwhite) ## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## HuntDay 19204.9 19204.9 1 188 11.0107 0.001088 **
## ID 3168.5 633.7 5 188 0.3633 0.873151
## HuntDay:ID 23341.2 4668.2 5 188 2.6764 0.023103 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Looking at all of the possible interactive effects, we can see that the only statistically significant effect is the interaction of hunting season day and how close the covey stayed in the shelter of hardwood forest.
summary(mixed_bobwhite)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: HW_Dist ~ (HuntDay * ID) + (1 | ID)
## Data: bobwhite
##
## REML criterion at convergence: 1969
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.85694 -0.39986 -0.09686 0.24302 3.07350
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 3679 60.65
## Residual 1744 41.76
## Number of obs: 200, groups: ID, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 16.6996 61.8382 188.0000 0.270 0.7874
## HuntDayRabbit 30.7117 14.0785 188.0000 2.181 0.0304 *
## ID2 -3.1075 87.2444 188.0000 -0.036 0.9716
## ID3 71.4500 87.7290 188.0000 0.814 0.4164
## ID4 -12.6545 87.8670 188.0000 -0.144 0.8856
## ID5 -0.6957 87.6185 188.0000 -0.008 0.9937
## ID6 -16.6996 87.8670 188.0000 -0.190 0.8495
## HuntDayRabbit:ID2 33.2485 19.1667 188.0000 1.735 0.0844 .
## HuntDayRabbit:ID3 -3.5830 21.1724 188.0000 -0.169 0.8658
## HuntDayRabbit:ID4 -31.1473 23.0762 188.0000 -1.350 0.1787
## HuntDayRabbit:ID5 -27.9975 22.5121 188.0000 -1.244 0.2152
## HuntDayRabbit:ID6 -22.5867 22.9181 188.0000 -0.986 0.3256
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) HntDyR ID2 ID3 ID4 ID5 ID6 HDR:ID2 HDR:ID3
## HuntDayRbbt -0.167
## ID2 -0.709 0.118
## ID3 -0.705 0.118 0.500
## ID4 -0.704 0.117 0.499 0.496
## ID5 -0.706 0.118 0.500 0.497 0.497
## ID6 -0.704 0.117 0.499 0.496 0.495 0.497
## HntDyRb:ID2 0.123 -0.735 -0.152 -0.086 -0.086 -0.087 -0.086
## HntDyRb:ID3 0.111 -0.665 -0.079 -0.183 -0.078 -0.078 -0.078 0.488
## HntDyRb:ID4 0.102 -0.610 -0.072 -0.072 -0.179 -0.072 -0.072 0.448 0.406
## HntDyRb:ID5 0.104 -0.625 -0.074 -0.074 -0.073 -0.162 -0.073 0.459 0.416
## HntDyRb:ID6 0.103 -0.614 -0.073 -0.072 -0.072 -0.072 -0.180 0.451 0.408
## HDR:ID4 HDR:ID5
## HuntDayRbbt
## ID2
## ID3
## ID4
## ID5
## ID6
## HntDyRb:ID2
## HntDyRb:ID3
## HntDyRb:ID4
## HntDyRb:ID5 0.382
## HntDyRb:ID6 0.375 0.384We can check our model using the Performance package:
performance::check_model(mixed_bobwhite)
There is an unequal number of points per day for each individual covey, which makes our study unbalanced. We can check the accuracy of this model and account for this by using emmeans:
bobwhite_means <- bobwhite %>%
group_by(HuntDay) %>%
summarise(mean_HW_Dist=mean(HW_Dist),
se_HW_Dist=sd(HW_Dist)/sqrt(n()))
bobwhite_means## # A tibble: 2 × 3
## HuntDay mean_HW_Dist se_HW_Dist
## <chr> <dbl> <dbl>
## 1 Quail 22.3 4.29
## 2 Rabbit 57.0 5.30bob_means <- emmeans(mixed_bobwhite, "HuntDay")## Cannot use mode = "kenward-roger" because *pbkrtest* package is not installed## NOTE: Results may be misleading due to involvement in interactionsbob_means## HuntDay emmean SE df lower.CL upper.CL
## Quail 23.1 25.3 188 -26.92 73.1
## Rabbit 45.1 25.1 188 -4.31 94.5
##
## Results are averaged over the levels of: ID
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95em_bob <- as.data.frame(bob_means)
em_bob## HuntDay emmean SE df lower.CL upper.CL
## 1 Quail 23.08173 25.34776 188 -26.920856 73.08431
## 2 Rabbit 45.11574 25.05572 188 -4.310749 94.54224em_bob$HuntDay <- factor(em_bob$HuntDay, levels= c("Quail","Rabbit"))
ggplot(em_bob, aes(x=HuntDay, y=emmean)) +
geom_point(size=5, color="#87b8d3") +
geom_errorbar(aes(ymin=emmean-SE, ymax=emmean+SE), width=.2, color="#6699CC") +
geom_point(size=5, data=bobwhite_means, x=bobwhite_means$HuntDay, y=bobwhite_means$mean_HW_Dist, color = "#336699") +
theme(axis.text.x = ggtext::element_markdown(color = "tan4", size = 12)) +
scale_x_discrete(labels = labels) +
theme(plot.caption=element_text(size=9, hjust=0, margin=margin(15,0,0,0)))+
theme_bw()+
labs(title="Comparing raw means to emmeans", caption="Light blue = raw means, Dark blue = adjusted means", x= "Hunting Day (1=Quail, 2=Rabbit)", y = "Average distance from hardwood cover (m)") 
We can see that the raw means and the adjusted emmeans don’t differ much, which lets us know that our model works for this analysis without having to account for the unbalanced points.
\(~\)
\(~\)
Nested Hierarchical Model
Great Tits (Parus major)
Advertisement signaling is usually linked to intersexual selection and intrasexual competition and thus is a key component of a species’ ecology. Using a novel spatial tracking system, the authors tested whether or not the spatial behavior of male and female great tits (Parus major) changes in relation to the response of a territorial male neighbor to an intruder. They tracked the spatial behavior of male and female great tits (N = 13), 1 hr before and 1 hr after simulating territory intrusions. The individual bird is a random effect, and the sex of the bird and measurements taken before and after are fixed effects.
tit <- read.csv("great_tits.csv")
tit$sex <- as.factor(tit$sex)
ggplot(tit, aes(sex, dist_m, colour = as.factor(id), shape=as.factor(measure))) +
geom_jitter(width =0.15, size=5, alpha=0.6)+
ylab ("Response Distance from Playback Location (meters)") +
xlab ("Sex") +
annotate("text", x = 2, y = 76, label = "13 birds") +
annotate("text", x = 2, y = 72, label = "2 measurements per bird") +
annotate("text", x = 2, y = 67, label = "26 total measurements", size=4, color="blue")+
labs(title="Sex-specific Responses to Territorial Intrusions", caption="1 = female
2 = male")+
scale_shape_discrete(name= "Measurements
Before / After")+
theme_bw()
In this study, 13 individuals were measured 2 times, 1 hour before playback and 1 hour after. The individuals are separated by sex to look at the different distances between females and males. Distance on the Y axis is measured in meters, with distance being how many meters away a bird was detected before and after a playback recording was played of a territorial male’s song. Based on the measurements of distance before and after playback, it looks like there isn’t really a pattern of whether or not sex plays a role in an individual bird’s movement after a territorial playback.
Using a linear mixed model to see if there’s an interaction between the sex of the bird and distance of detection:
lmmtit <- lmer(dist_m ~ (1|sex/id), data = tit)
summary(lmmtit)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dist_m ~ (1 | sex/id)
## Data: tit
##
## REML criterion at convergence: 226.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3975 -0.5381 -0.1098 0.4360 1.6636
##
## Random effects:
## Groups Name Variance Std.Dev.
## id:sex (Intercept) 236.40 15.38
## sex (Intercept) 0.00 0.00
## Residual 99.41 9.97
## Number of obs: 28, groups: id:sex, 14; sex, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 27.711 4.521 13.000 6.13 3.6e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularanova(lm(dist_m ~ sex/id, data = tit))## Analysis of Variance Table
##
## Response: dist_m
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 47.4 47.40 0.1309 0.7207
## sex:id 2 90.6 45.31 0.1251 0.8830
## Residuals 24 8692.3 362.18Our p-values are not reading as statistically significant, which matches what we can see in the plot. What this can be interpreted as is that there is not a significant interaction between sex of the bird and the distance it can be detected before or after a territorial playback is used.
